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Abstract: Inaccurate description of the equilibrium can yield to spurious effects in gyrokinetic 
turbulence simulations. Also, the Vlasov solver and time integration schemes impact the conser- 
vation of physical quantities, especially in long-term simulations. Equilibrium and Vlasov solver 
have to be tuned in order to preserve constant states (equilibrium) and to provide good conser- 
vation property along time (mass to begin with). Several illustrative simple test cases are given 
to show typical spurious effects that one can observes for poor settings. We explain why Forward 
Semi-Lagrangian scheme bring us some benefits. Some toroidal and cylindrical GYSELA runs are 
shown that use FSL. 
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Precision du mouvement de particules non perturbees dans 

un code gyrocinetique 

Resume : Une description imparfaite de l'equilibre peut conduire a des artefacts numeriques 
dans des simulations turbulentes gyrocinetiques. Aussi, le solveur de Vlasov and le schema 
d'integration en temps ont un impact fort sur la conservation de quantites physiques, notamment 
lorsque les simulations sont en temps long. L'equilibre et le solveur Vlasov doivent etre finement 
choisis et parametres pour conserver les etats constants (equilibre) et pour autoriser de bonnes 
proprietes de conservation en temps (en commengant par la conservation de la masse). Plusieurs 
cas tests illustratifs sont donnes pour montrer les problemes numeriques typiques que Ton peut 
observer si les choix pris ne sont pas adequats. Nous expliquons pourquoi le schema FSL 
(Forward Semi-Lagrangian) apporte une reponse a certains problemes. Des simulations en 
configuration cylindrique et torique de GYSELA sont presentees qui utilisent FSL. 

Mots-cles : gyrocinetique, lois de conservation, FSL 



Unperturbed motion of particles 



3 



1 Introduction 

Inaccurate description of the gyrokinetic equilibriu m can yield unphysical excitation of zonal 
flow oscillations [ABH + 06]. Moreover, as stated in |IIK + Q8] . it is important to define the initial 
condition using a relevant gyrokinetic equilibrium, especially in collisionless full-f simulations. 
In the following, accuracy aspects are investigated for both the gyrokinetic equilibrium and 
the Vlasov solver used in the GYSELA code. If proper care is not taken for Vlasov solver 
and gyrokinetic initial equilibrium, one can observe that some conservation properties are not 
satisfied, for example mass conservation. 

The gyrokinetic framework of the study is explained in the next section. Then, quite exper- 
imental investigations are shown to understand numerics associated with equilibrium and mass 
conservation. Forward Semi-Lagrangian (FSL) scheme is presented briefly. The splitting in a 
linear and a non-linear part of the Vlasov equation is explained. Results of numerical simulation 
for 4D cylindrical test cases and toroidal test cases are presented in the last section. Forward 
Semi-Lagrangian scheme helps a lot to preserve good accuracy on mass and total energy. 



2 Description of the context 
2.1 Gyrokinetic Vlasov equation 

The time evolution of the gyro-center coordinates (x, v\\ , ji) of species s is given by the collision- 
less electrostatic gyrokinetic equations: 

dx 1 - 

— = v\\b* a • Vx l + v Esgc • Vx l + v Ds • Vx l (1) 



II ^ -> - m s v\\ -+ , x 

ms dt = ' VB " 6 A ' W + B * ESGC ' VB (2) 



where x % corresponds to the z-th covariant coordinate of x, B is the magnetic field and B^ s and 
6* are defined as: 

m s v\\ -> -> 

B; s = B+^-^ b-J (3) 
S Bis e.Bl, B [ } 

The advection terms are: 

t* . ^ _ u» _ B-Vx^ m a v n Mo J ■ W 

s ~ s ~ + B (5) 

\\s S \\s 

v Ds • Vx* = v % Da = i e ^ B J [5, x*] (6) 

vesoc • Vx* = v% SGC = — [0, x l ] (7) 

^11- 

0EaGc'VB = --^r[B,4>] (8) 

II -s 

the Poisson bracket is defined by [f,g] = 6-(vfxVG). The velocity parallel to the magnetic field 
is vn. The magnetic moment \i = m s v\/ (2B) is an adiabatic invariant with v_\_ the velocity in 
the plane orthogonal to the magnetic field. Vacuum permittivity is denoted /zq. The term ve gc 
represents the electric E x B drift velocity of the guiding-centers and ve> the curvature drift 
velocity. The Jacobian in phase space is J 3 (r,Q)B*\ 3 (r,0,v\\). 

Some references concerning the framework we use to solve this equation can be found in 
jGBB+061 lGS?F08l IBCG+Ill I SGA+Ill IAGG+11) . 



2.2 Quasineutrality equation 

The quasineutrality equation and parallel Ampere's law close the self-consistent gyrokinetic 
Vlasov- Maxwell system. However, in an electrostatic code, the field solver reduces to the 
numerical solving of a Poisson-like equation [Hah88j . In tokamak configurations, the plasma 
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quasineutrality (denoted QN) approximation is currently assumed [GBB+06]. Electron inertia 
is ignored, which means that an adiabatic response of electrons is supposed. We define the 
operator Vj_ = (d r ,^do). We note no the equilibrium density, Bo the magnetic field at the 
magnetic axis and T e (r) the electronic temperature. We have also B(r,0) the magnetic field, 
Jo the Bessel function of first order and k± the transverse component of the wave vector. 
Hence, the QN equation can be written in dimensionless variables 



n (r) 

where p is defined by 



n (r) 



V±$(r,0,<p) 



T e (r) 



[Q(r,0,<p)-(*) FS (r)] = p(r,0,<P) 



p(r, 0, ip) = J B^ s dp J dv\\J (k ± ^2p)(f - f eq )(r,6,ip,v h p). 



(9) 



(10) 



with f eq representing local ion Maxwellian equilibrium, and (.) FS (r) the average on the flux 
surface labelled by r. 



3 Unperturbed motion of particles 
3.1 Simplification of the equations 

Let suppose that a collisionless equilibrium distribution function f eq is taken for the initial dis- 
tribution function: f t =o(r, y?, v\\). By definition, if the Vlasov and Quasi- neutrality equations 
are accurately solved, the distribution function f t should remain constant over time, i.e. equal 
to f eq . Also, the electric potential <£(r, #, tp) should remain equal to zero because the right hand 
side of the Quasi- neutrality equation, the integral of / — f eq , is zero. 

Then, one quantitative measure of Vlasov solver accuracy is its capability to preserve an equilib- 
rium distribution function over time, assuming $ = 0. In order to get a simple tool to evaluate 
the Vlasov solver, one can also restrict to the case \i — 0. It is worthwhile rewriting gyrokinetic 
equations in this reduced setting: 

dr i . . 

— =v\\b* a .Vx t + v Da .Vx % (11) 

d^ii , N 

m s -^=0 (12) 



dt 



It gives 



dx l B • Vx % m s v\f n J • Vx 
-v\ 1 



(14) 



dt B* s 11 e s B* s B \e s B* {s B 
dv\ 
"dt 

and finally, assuming that the current J is perpendicular to the poloidal plane 

„2 



[B,j] (13) 



dr B r m 8 vt 



-V\\ 



[B,r] (15) 



dt 5* 11 e s B* B 

\\s s \\s 

\\s S \\s 

dtp _B m s «|| 2 J v m s vj 

^-Bf;* + ^l3 + ^B^ ^ 

Particles keep a constant parallel velocity over time. Then the trajectory of one particle in 
phase space remains in the 3D space dimensions. 



3.2 Trajectory of particles 

The GYSELA code uses a Semi-Lagrangian scheme combined with a Strang splitting. 
To guarantee second order accuracy in time, the following sequence [GBB+06] is used: 
(v\\/2,(p/2,r0,(p/2,v\\/2), where factor 1/2 denotes a shift over At/2. In this sequence, v\\ 
and shifts correspond to ID advections along v\\ and tp directions, and rO represents a 2D 
advection (displacement in r and are strongly coupled). Let us now focus on the discretization 
of the 2D advection which requires specific attention. 
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Taylor expansion For a single rO advection in the semi-Lagrangian approach, a displace- 
ment (Sr,S0) is computed to find the foot of the characteristics that ends at a grid node 
(ri,9j). This disp lacement is computed thanks to a Taylor expansion of the fields acting on 
particles [GBB+06]. One can derive two fields a and /3 from Eqs. 

(ll) = a(n, 3 )At + (3(r u 3 )At 2 + 0(At 3 ) 

The problem is that, if Sr or SO is larger than grid element sizes denoted Ar and A0, then this 
approximation can be really inaccurate. 



15| 16| 17 ) such as: 



Precomputed trajectories Another method is described here to find out the foot of the 
characteristics associated to each grid node. Since the fields acting on particles, Eqs. ( 11|12 ), 
do not depends on time t, one can approximate (Sr, S0)^ r .^.^ once for all. We have used Runge- 
Kutta time integration scheme (RK2) with a small time step to precompute these trajectories. 
This approach is possible only for the linear terms of Eqs. ( 11|12 ), but not for non-linear terms, 
present in Eqs. (T|2), that depend on E t and Let us choose a St, such as M St = At with 
MgN and M large enough. One can build a series r n and n as follows: 



St 



a(r n ,0 n ) 



(^i) = (^)+««(^ i > » n+i ) 

The initial condition is set to r° = r^, 0° = 0j. After M steps of Runge-Kutta iteration, it gives 




Because these trajectories can be computed only once when the simulation starts, M can be 
taken quite large (we will assume M = 64 in the following). These precomputations will not 
impact significantly the global simulation time. Other time integration schemes have been tried: 
RK3, RK4, and also larger values of M, but these modifications does not impact significantly 
the precision in our test cases. 



3.2.1 Forward Semi-Lagrangian scheme 

The Forward Semi-Lagrangian scheme [CRS09J (the acronym is FSL) is an alternative to the 
classical Backward Semi-Lagrangian scheme (BSL). It is based on forward integration of the 
characteristics. The distribution function is updated on an eulerian grid, and the pseudo- 
particles located on the mesh nodes follow the characteristics of the equation forward for one 
time step, and are then deposited on nearest nodes. While the main cost of the BSL method 
comes from the interpolation step, FSL spends most of computational time in the deposition 
step. The FSL scheme can be set up such as mass is preserved along time. 

We present the approach taken for the 2D advection in (r, 0) . In our setting, the algorithm 
of Forward Semi-Lagrangian is the following (with S the cubic B-spline): 

• Step 1: Compute the spline coefficients l (f n is known) such that 
(J s Bl la f n ) id = J2<iS(xi - x k )S( yj - Vl ) 

k,l 



• Step 2: Integrate forward in time the characteristics from t n to t n+1 , given as initial data 
the grid points (xk,yi)- We obtain the new particle positions (x%,yf) at time t n+1 by 
following the characteristics. 

• Step 3: Deposition of particles localized at (x^,yf) 

fit 1 = ( E w h s ^ - x k) s (y 3 - vl) ) /(J*B\\a)ij 
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3.3 Numerical experiments 

3.3.1 Accuracy of conservation of invariants 

An equilibrium solution of the collisionless gyrokinetic equation must satisfy some conditions. 
To get an equilibrium for the Vlasov equation, it suffices to take an arbitrary function of con- 
stants of motion in the unperturbed characteristics. In an axisymmetric toroidal configuration, 
a gyrokinetic Vlasov equilibrium is defined by three constants of motion: the magnetic mo- 
ment p, the energy ^vJ/2+^B(f,^), and the canonical toroidal angular momentum |AGG + ll) 
p v =<0(r)+/v||/.B(r,0). The ip(r) function is defined thanks to the safety factor q(r) by the relation: 

dip I dr— — Bo r/q{r). 

To get a good accuracy in integrating the Vlasov equation in time, and then to conserve 
exactly an initial equilibrium distribution function, we have to take care of the correctness of: 

• the initial equilibrium setting; distribution function will be a function of the constants of 
motion: p,£,P<^; 

• computing the displacements in the Vlasov solver (finding the foot of the characteristic); 
a sufficiently small time step has to be chosen for the time integration scheme for example, 
but also the resolution in space must be sufficient to have enough accuracy to interpolate 
displacement fields; 

• interpolate the distribution function in the Vlasov solver with enough resolution; the 
number of points in the computational domain has to be tuned for this issue; 

We will be able to estimate quantitatively the quality of the distribution function after several 
time steps by: 

• measuring the difference between the initial equilibrium function and the distribution 
function at a given time step. 

3.3.2 First test case: equilibrium conservation 

The initial distribution function is set up to 

f t=0 (r,e,<p,v ll ) = l [PvuPv2] (P v -P vl ) m (P v -P v2 rsin( 7 P ip ) (18) 
where the invariant P^ is a function of r, vu and 1 is the indicator function. We will look at 
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Figure 1: Initial function f t =o at v\\=—3vth0i<p = (left), v\\=0,p = (right) 
the norms u p (t) = \\ft — /t=o||p with p = 1, oo. They are defined as 

iZoo(t) = sup{\f t - f t =o\ : V(r,0,<p,V||)} 
ui(t) = J \f t - f t =o\J s B*\ s drdOd(pdv\\ 

Let us take the following parameters : 7 = 47r/(P^,2 — P<pi), ra = 1, p* = 0.01. The values 
(P^i, P<p2) are choosen adequately in order to have an interval [P^i, P^} that is strictly included 
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in the computational domain r G [r m i n ,r max ],0 G [0, 27r], i>|| =-3%o- A poloidal cut of this 
initial function is sketched for two values of parallel velocity in Figure [T] As the initial function 
depends only on the invariant P^, it should be conserved by time integration of the Vlasov 
equation. The duration of the experiment is denoted t max . The timing unit is the ion cyclotronic 
time f2 _1 , it will also be the case in all plots of this document. The Vlasov equation is solved 
during all time steps with forced to zero. An ideal simulation should preserve all norms 
Vp, u p (t) = 0. Practically and in the worst case, the function u p (t) will move quickly away 
from as t grows. 
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Figure 2: Evolution of the L\ norm u\ (t) with Figure 3: Evolution of the norm Uoo(t) 
N r = 512, N e = 512 with N r = 512, N e = 512 
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Figure 4: distance between f t and f t =o (modulus of the difference) at t = 100 with At = 10 
and Precomputed trajectories (left), distance \\f t — ft=o\\ at £ = 100 with At = 10 and Taylor 
expansion (right) 



In the experiment of Figures [2] and [3j the Li, norms are drawn for some simulations. A 
comparison is made between the two methods presented previously that are used to calculate 
the feet of the characteristics. The parameter At is set to 2 and 10, which is near typical values 
used in GYSELA simulations and t max = 100. The computational domain is quite refined in 
(r, 6) with N r = 512, Nq = 512. One can remark that Precomputed trajectories strategy is 
more accurate than Taylor expansion. Nevertheless, the difference shown here is significant for 
At = 10, but quite negligible for small At = 2. On Figure [4j the differences to the initial 
distribution function at final state, are shown for an arbitrary value of v\\ = —Sv t ho- One can 
see here where errors are located using both strategies: Precomputed trajectories at the center, 
and Taylor expansion at right. One can see that the amplitude of the error is larger in the case 
of Taylor Expansion, and also less spatially localized. 
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Figure 5: Evolution of the L\ norm ui(t) for Figure 6: Evolution of the norm ^oo(t) for 
different poloidal mesh sizes different poloidal mesh sizes 



On Figures [5] and [6) some L-p norms produced by three simulations are shown (Precomputed 
trajectories method is used). Between two consecutive simulations, the sizes of the domain 
(N r ,No) have been multiplied by a factor two. While increasing (N r ,Ne), one can notice that 
simulations become more and more accurate, i.e. L-p norms are reduced. The increase rate of 
the error depends on the factor N r Nq. 
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Figure 7: Evolution of the L\ norm u\(t) for Figure 8: Evolution of the 
different poloidal mesh sizes different poloidal mesh sizes 



,(«) for 



The domain sizes A^ r and Nq are increased independantly each other in the Figures [7] and 
[8] (Precomputed trajectories method is used). It is clear that the error is only of the order of 
the smallest size among N r and Nq sizes. As far as Equations ( 1 1|12 ) are concerned, it is then 
noteworthy that 7V r and Nq should be chosen close to each other in order to improve accuracy. 



3.3.3 Second test case: shear flows 

Compared to the previous experiments, the initial distribution function considered in the present 
section is taken with an inhomogeneity in 0. Only the Precomputed trajectories method has 
been used for these simulations, because it is more accurate in any case. As variable 6 is not 
an invariant, the distribution function is expected to evolve in time. For the sake of simplicity 
and to reduce computational costs, the simulation is setup with N^ = l,iV V|| = 1. The initial 
distribution function is now 



](P<p 
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Figure 9 
/ at t = 
(right) 



: Initial function f t =o with vn = —3 v t ho, N r — 512, Nq = 512 (left), distribution function 
2000 with BSL scheme (middle), distribution function / at t = 50000 with BSL scheme 
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Figure 10: Initial function f t= o with vn = —3vtho,N r = 128, iV^ = 128 (left), distribution 
function / at t = 2000 with BSL scheme (middle), distribution function / with BSL scheme at 
t = 50000 (right) 
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Figure 11: Initial function f t= o with vn = —3vtho,N r — 128, iV^ = 128 (left), distribution 
function / at t = 2000 with FSL scheme (middle), distribution function / at t = 50000 with 
FSL scheme (right) 



On Figures [9| [lO 11, the distribution function is shown for different time steps (parallel 
velocity is v\\ = — 3 vtho, ra is set to 1, , p* = .01). One can notice the shearing that appears along 
time. The system dynamics tends to create filamentation in configuration space. Nevertheless, 
the mass ||/||i = J f J s Bidrd9d</)dv\\ and the maximum value ||/||<x> = su P{\f\ V(r, 0, 0, v\\)} 
should theoretically be conserved at any time. Finer structures than mesh element size are 
dynamically created (see the rightmost plot of Figure [9]), and they can not be correctly captured 
(filamentation induces that mass can be lost or gained at each time step). As it can be seen in 
Figures 12 and 13 , with the classical Backward Semi-Lagrangian scheme ||/||oo quantity does not 
remain constant, ||/||i is well conserved at (iV r = 512, N e = 512) while (iV r = 128, 7V# = 128) 
does not lead to good conservation of ||/||i. It is noteworthy, that finer resolution in space - 
going from (N r = 128, iV^ = 128) to (7V r = 512, iV^ = 512) greatly improves the conservation 
of mass and of the largest absolute value. If one takes the Forward Semi-Lagrangian scheme, 
the Figures 12 and 13 illustrate that ||/||oo behaves like BSL but ||/||i is better preserved, as 



expected. Then, FSL scheme is a good solution as far as mass conservation is required. 
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Normalized mass of f (Nr=128 Ntheta=128) - FSL 
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Figure 12: Time evolution of the mass of / Figure 13: Evolution of the / maximum value 
for different poloidal mesh sizes for different poloidal mesh sizes 



Let suppose one wants to keep the framework of the Backward Semi-Lagrangian scheme 
and to avoid FSL, another solution to preserve mass could be to have an aligned mesh. If 
one change the variables (r, 9,p : v\\) to (P^, 0, tp, v\\) : the two initial distribution functions we 
have seen are more accurately discretized. Also, the conservation of mass is improved for the 
following reason: advection is done only in the direction there is no shift along the P^ 
direction. The displacement in the poloidal plane during one time step amounts to shift only in 
0. So mass has just to be conserved on each mesh line labelled by P^, and it will automatically 
give global mass conservation. It is easy to construct an advection operator and the associated 
interpolation scheme that preserve mass on each mesh line P^. Another benefit is that in 
this new set of variables, the equilibrium function is better discretized and represented: steep 
gradient is only in P^ direction in the poloidal plane. Nevertheless, building and managing a 
mesh that includes P^ variable is much more difficult to implement. 



3.3.4 Extrapolation to nonlinear test cases 



The first test case of section |3.3.2| is linked to what is going on at the beginning of a Gysela 
run. The initial distribution function is a perturbed equilibrium. One wishes to remain close 
to this quasi-equilibrium (if the amplitude of electric potential is very small) without beeing 
too much biased by some large numerical artifacts of the linear operator that we have seen: 
artificial increase of ui : i^. 



The second test case of section 3.3.3 is near what is happening in non-linear dynamics in 
which many localized spatial structures develop and grow. One expect to conserve the total 
particle density, even if the distribution function becomes something else than an ideal function 



of motion invariants. The decrease of maximal value of Figure 13 is not acceptable to perform 
accurate simulations. Aligned coordinates along P^ should be able to improve conservation of 



4 Nonlinear gyrokinetic simulations 



4.1 Split linear and nonlinear parts 

4.1.1 Global separation of linear /nonlinear terms 



The equations (TJ2) can be split into two parts (in the same spirit that |IIK + Q8p . The first 
part includes the nonlinear terms that depends on the electric potential. The second part 
comprises the other terms. One can solve these two parts one after the other by splitting. 
On the first hand, the linear operator is presented in Eqs ( 19|20 ), on the second hand, the 
nonlinear operator is described by Eqs ( 21|22 ). 
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Linear operator C 



Nonlinear operator J\f 



dx l 



dt 

dv\\ -> -> 

m s = — [ib* • VB 

dt 



vnb* • Vx l + v Ds • Vx x 



(19) 
(20) 



dx l 
~dt 



- V Es 



■ Vx l 



m s v\\ _ 



(21) 
(22) 



The linear operator exhibits large displacements at large modulus of parallel velocity, and 
also induces shear flows. These features can interfere with the nonlinear dynamics that even- 
tually involves small displacements that are not of the same order of magnitude. Moreover, 
as the dynamics generated by the two operators are different, induced accuracy problems have 
possibly not the same characteristics for the two operators; and the limitations (CFL-like con- 
ditions) on the time step are also not the same. Ideally, one should be able to fix the time step 
of linear and nonlinear operators independantly in order to enhance accuracy. 

The original formulation that uses Strang splitting in phase space for the Vlasov solver, 
but without a linear /nonlinear splitting, can leed to serious troubles at high parallel velocities. 
Indeed, directional splitting suffers from some large shifts in ip, r and 6 at each directional ad- 



vection (section 3.2). Then, for not so large At, the evaluation of electric fields \ is not done 
at the right spatial position at each substep of the splitting except for the first substep of the 
splitting. One has to take very small time step to recover small displacements in the linear 
operator and then a reasonable accuracy in the evaluation of E^y 

The splitting between linear and nonlinear parts eliminates this problem. The nonlinear oper- 
ator is applied alone, thus the linear operator and its large shifts does not interact badly with 
nonlinear solver. The proposed solution is to perform at each time step: first, a directional split- 
ting for the linear operator C with the sequence (Vy/2, <£/2, r#, <£/2, v\\/2), second, a directional 
splitting for the nonlinear operator M with the same sequence (v\\/2,(p/2,rQ,(p/2,v\\/2). 

This approach is a liitle bit expensive in term of computational cost because it nearly doubles 
time to solution. 



4.1.2 Conservation of mass 

Each of the operators £, J\f conserves the mass independently. One can track during simulation 
which operator degrades or not the total mass. 

4.2 Experiments 

In the following GYSELA runs, the BSL and FSL schemes are compared for two test cases 
with \i =0 (4D only). Results for both BSL scheme and FSL scheme are shown. In the FSL 
version that has been run, FSL scheme is present in all advections: ID (ty?,^) and 2D (r, 6). 
The rest of the code is unchanged compared to the BSL version. In order to simplify boundary 
conditions in v\\, the distribution function is fixed at v m i n and v max . Therefore, we have also 
taken this choice for the BSL runs presented hereafter. 



4.2.1 Cylindrical case 
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Figure 14: Evolution of the mass for a 4D, Figure 15: Evolution of the energies for a 
\i = 0, cylindrical test case 4D, \i = 0, cylindrical test case 
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Some plots for our reference test case data4D_slab_A32 (p* = 1/32, p = 0) is given in 
Fig. ( |14|15|21|22| . The computational domain size is N r = 128, N e = 256, = 32, A^, = 48. 

The mass is well preserved with FSL as it can be seen on Fig. [l4j Potential energy, kinetic 
energy and total energy are presented on FigfTK) The total energy should theoretically remains 
at zero, the two schemes are quite proximate to that. During the beginning of the simulation, 
FSL and BSL give quite similar electric potential as it is shown at time t = 1920 in Fig. ( |21|22 ). 
Nevertheless, after a while, the two schemes diverges slowly (see for example t = 8000). 
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Figure 16: Poloidal cut of electric potential at different time steps (t = 1920 left, t = 8000 
right), (f = 0, BSL scheme 
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Figure 17: Poloidal cut of electric potential at different time steps (t = 1920 left, t = 8000 
right), <p = 0, FSL scheme 



4.2.2 Toroidal case 



A small toroidal test case has been set up (p* = 1/40). The computational domain size is 
N r = 256, N e = 256, = 64, A^, = 48. The FSL version outperforms the BSL version from 
the point of view of mass conservation. The conservation of total energy is far more better with 
FSL and this is due to improved mass conservation. A precise tracking of mass losses/gains 
allows us to conclude that the linear operator C is mainly responsible for mass degradation 
using BSL approach. An explanation is that in toroidal setting, C causes large displacements, 
whereas in cylindrical setting (previous subsection) this was not the case. 
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4D toroidal case - mu=0 
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Figure 18: Evolution of the mass for a 4D, \i = 0, toroidal test case 
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Figure 19: Evolution of the energies for a 
4D, [i = 0, toroidal test case with FSL 
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Figure 21: Poloidal cut of electric potential 
right), ip = 0, BSL scheme 
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Figure 20: Evolution of the energies for a 
4D, [i = 0, toroidal test case with BSL 
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Figure 22: Poloidal cut of electric potential at different time steps (t = 
right), cp = 0, FSL scheme 
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5 Conclusion 

Few solutions have been proposed to help preserving constant states, L\ and ^ norms in 
gyrokinetic semi-Lagrangian simulations. Some GYSELA runs show that we get benefits from 
these improvements in term of accuracy. 
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